%% Table4.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Date: May 4, 2019
%
% Program: this is a first attempt at operationalizing the granular
% residual calculations following from master_notesv3.pdf (page 34) using
% the methodology of Gabaix and Koijen's "Granular IV" paper. To begin, we
% will do this based on PWT data for real GDP growth as well as TFP
% measures reported in this dataset. 
%
% The PWT dataset will be for 1970-2014, and will be loaded as levels, and
% then manipulated.
%
% Once the granular residuals are constructed based on the different
% partner country growth rates, we will compare to France's actual real GDP
% and TFP growth to see if any patterns exist. I.e., the model generated
% elasticities that can map the country shocks into an aggregate foreign
% component that can predct overall French GDP or TFP growth.
%
% We can also compare to the granular residuals we've calculated from
% previous work using French firm-level data, but will need to restrict the
% sample size over years further for this.
%
% UPDATED 2/21/20 (JDG): adjusted calls of directories given still in
% Robustness folder. Will need to be adjustd later. Also switched to double
% deflation method for calculating elasticities.
%
% Also created a version of Table 4 of paper to output.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%
% Loan PWT real gdp, exchange rate, and TFP data from 1970-2014

% rgdpe   Expenditure-side real GDP at chained PPPs (in mil. 2011US$)
% rgdpo   Output-side real GDP at chained PPPs (in mil. 2011US$)
% rgdpna  Real GDP at constant 2011 national prices (in mil. 2011US$)
% ctfp    TFP level at current PPPs (USA=1)
% cwtfp   Welfare-relevant TFP levels at current PPPs (USA=1)
% xr      Exchange rate, national currency/USD (market+estimated)
% 
 load MAT/External/pwt9_rgdpftp.mat
% 
 ISO = string(pwt9rgdptfp.ISO);
 year_pwt = pwt9rgdptfp.year;
% rgdpe = pwt9rgdptfp.rgdpe;
% rgdpo = pwt9rgdptfp.rgdpo;
% rgdpna = pwt9rgdptfp.rgdpna;
 ctfp = pwt9rgdptfp.ctfp;
 cwtfp = pwt9rgdptfp.cwtfp;
% xr = pwt9rgdptfp.xr;


load MAT/External/RGDPlcu_WDI.mat

ISO_WDI = string(RGDPlcuWDIclean.ISO);
year_WDI = RGDPlcuWDIclean.year;
rgdpe = RGDPlcuWDIclean.Value;


%%
% Caculate growth rates of rgdpe (in local currency_ and ctfp
% Will fill in NaN for missing values and keep a balanced panel for now

% Transform RGDP to local currency to get rid of exchange rate swings for

%rgdpe_lc = rgdpe.*xr;
rgdpe_lc = rgdpe;
% A numerical index to loop over countries

id = grp2idx(ISO);
id_WDI = grp2idx(ISO_WDI);

france = 15;

% Initialize vector of growth rates for real GDP and TFP

gr_rgdpe_lc = NaN(size(id_WDI,1),1);
gr_ctfp = NaN(size(id,1),1);

for i = 1:max(id)
   
   x = ctfp(id==i);
   grx = NaN(size(x,1),1);
   grx(2:end) = x(2:end)./x(1:end-1)-1;
   gr_ctfp(id==i) = grx;   
   
   clear x grx;
end

for i = 1:max(id_WDI)
   x = rgdpe_lc(id_WDI==i);
   grx = NaN(size(x,1),1);
   grx(2:end) = x(2:end)./x(1:end-1)-1;
   gr_rgdpe_lc(id_WDI==i) = grx;
   
   clear x grx;
   
end

%%
% Calculate Granular residual based on partner countries:
%
%  a) TFP growth
%  b) Real GDP growth
%
% Note: the numnber of observations may vary per country, and we
% will compare to overall GDP and TFP growth of France later.

CTY = {'AUS', 'AUT', 'BEL', 'BGR', 'BRA', 'CAN', 'CHN' ,'CYP', ...
'CZE', 'DEU', 'DNK', 'ESP' , 'EST', 'FIN', 'FRA', 'GBR', ...
'GRC', 'HUN', 'IDN', 'IND', 'IRL', 'ITA', 'JPN', 'KOR', ...
'LTU' ,'LVA', 'MEX','MLT', 'NLD', 'POL', 'PRT', 'ROU','RUS', ...
'SVK', 'SVN', 'SWE', 'TUR', 'TWN', 'USA','ROW'};
% Initialize vectors we will collect results in

T = size(ctfp,1)/max(id);
T_WDI = size(rgdpe_lc,1)/max(id_WDI);
Gamma_nRGDP = zeros(T_WDI,size(CTY,2));
Gamma_nTFP = zeros(T,size(CTY,2));
Overall_nRGDP = zeros(T_WDI,size(CTY,2));
Overall_nTFP = zeros(T,size(CTY,2));

% Run loop a calculate Granular Residula w.r.t. to each country's shock. 
% NOTE: only do w.r.t. TFP growth now given we need to adjust formula for
% real GDP growth calculation -- imperfect what we do for GDP

for i = 1:size(CTY,2)
   cty = string(CTY(i));
   if cty~='FRA' & cty~='ROW'
   % Get required counterfactual data
   
   fname = strcat('MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_5_output_alpha4_rho3_',cty,'prod_p10_approx.mat');   
   load(fname,'VA_fj_0','VA_hat_fj','RGDP_cpi_hat_n','RGDP_def_hat_n','GDP_deflator_n','DoubleDeflation_deflator_n','RGDP_doubledef_hat_n');   
   
   % Calculate shares for weighted and unweighted aggregation
   
   w_sh = VA_fj_0/sum(VA_fj_0);
   uw_sh = 1/size(VA_fj_0,1);
   
   % Calculate elastitcity at firm level given 10% productivity shock to
   % country in counterfactual
   
   e_f = (VA_hat_fj-1)/.1; % F*1
   e_f_real = (VA_hat_fj/DoubleDeflation_deflator_n(france)-1)/.1;          % F*1
   e_f_GDP = (VA_hat_fj/DoubleDeflation_deflator_n(france)-1)./(RGDP_doubledef_hat_n(i,1)-1); 
%    e_f_real = (VA_hat_fj/GDP_deflator_n(france)-1)/.1;          % F*1
%    e_f_GDP = (VA_hat_fj/GDP_deflator_n(france)-1)./(RGDP_def_hat_n(i,1)-1);
   e_f(isnan(e_f))=0;
   e_f_real(isnan(e_f_real))=0;
   e_f_GDP(isnan(e_f_GDP))=0;
    cty
   mean(e_f)
   mean(e_f_real)
   mean(e_f_GDP)
   % Calculate Granular Residual component per country per time period. 
   % Given we are summing over firms and countries, we will sum over firms
   % per country in loop and then we will sum across countries after loop
   % to get total Granular Residual

   % (a) TFP Growth Granular Residual
   
   growth =  gr_ctfp(ISO==cty);      % T*1
   
   for t = 1:T
      granres_n = sum(w_sh.*e_f_real*growth(t)) - mean(e_f_real*growth(t));
      Gamma_nTFP(t,i) = granres_n;
      Overall_nTFP(t,i) = sum(w_sh.*e_f_real*growth(t));
      Mean_nTFP(t,i) = mean(e_f_real*growth(t));
      Mean_nTFP_(t,i) = Overall_nTFP(t,i)-Gamma_nTFP(t,i);
      clear granres_n;
   end
   
   clear growth
   
   % (b) Real GDP Growth Granular Residual
   if cty~='TWN'
   growth =  gr_rgdpe_lc(ISO_WDI==cty);      % T*1
   
   for t = 1:T_WDI
      granres_n = sum(w_sh.*e_f_GDP*growth(t)) - mean(e_f_GDP*growth(t));
      Gamma_nRGDP(t,i) = granres_n;
      Overall_nRGDP(t,i) = sum(w_sh.*e_f_GDP*growth(t));
      Mean_nRGDP(t,i) = mean(e_f_GDP*growth(t));
      Mean_nRGDP_(t,i) = Overall_nRGDP(t,i)-Gamma_nRGDP(t,i);
      clear granres_n;
   end
   end
   clear cty fname w_sh uw_sh e_f growth
   end
end

year_WDI_=unique(year_WDI);
year_pwt_=unique(year_pwt);

Gamma_TFP = nansum(Gamma_nTFP,2);
Gamma_RGDP = nansum(Gamma_nRGDP,2);
Overall_TFP = nansum(Overall_nTFP,2);
Overall_RGDP = nansum(Overall_nRGDP,2);
Mean_TFP = nansum(Mean_nTFP,2);
Mean_RGDP = nansum(Mean_nRGDP,2);
Mean_TFP_ = nansum(Mean_nTFP_,2);
Mean_RGDP_ = nansum(Mean_nRGDP_,2);

std(Overall_TFP(2:end))
std(Mean_TFP(2:end))
std(Mean_TFP_(2:end))
std(Gamma_TFP(2:end))

grTFP_FRA = gr_ctfp(ISO=="FRA");
grRGDP_FRA = gr_rgdpe_lc(ISO_WDI=="FRA");

% Summary statistics and figures
ratio_Gamma_TFP=Gamma_TFP(12:end)./grRGDP_FRA(2:55);
ratio_Gamma_RGDP=Gamma_RGDP./grRGDP_FRA;
nanmedian(ratio_Gamma_TFP)
nanmedian(ratio_Gamma_RGDP)
nanmean(ratio_Gamma_TFP)
nanmean(ratio_Gamma_RGDP)


b = regstats(Gamma_TFP(12:end),grRGDP_FRA(2:55),'linear','yhat');
figure(1)
plot(grRGDP_FRA(2:55),Gamma_TFP(12:end),'bo',grRGDP_FRA(2:55),b.yhat,'r-','LineWidth',2)
xlabel('GDP Growth')
ylabel('TFP Granular Residual')
saveas(gcf,'Figures/rho3eta1.000001lambda1.000001_psi3_sigma1.5/TFPGran_vs_GDPGrowth.pdf')
clear b
std(Gamma_TFP(12:end))
std(grRGDP_FRA(2:55))
corr(Gamma_TFP(12:end),grRGDP_FRA(2:55))
std(Gamma_TFP(42:58))
std(grRGDP_FRA(32:48))
corr(Gamma_TFP(42:58),grRGDP_FRA(32:48))

b = regstats(Gamma_RGDP(2:end-1),grRGDP_FRA(2:end-1),'linear','yhat');
figure(2)
plot(grRGDP_FRA(2:end-1),Gamma_RGDP(2:end-1),'bo',grRGDP_FRA(2:end-1),b.yhat,'r-','LineWidth',2)
xlabel('GDP Growth')
ylabel('GDP Granular Residual')
saveas(gcf,'Figures/rho3eta1.000001lambda1.000001_psi3_sigma1.5/GDPGran_vs_GDPGrowth.pdf')
clear b
std(Gamma_RGDP(2:end-1))
std(grRGDP_FRA(2:end-1))
corr(Gamma_RGDP(2:end-1),grRGDP_FRA(2:end-1))
std(Gamma_RGDP(15:end-1))
std(grRGDP_FRA(15:end-1))
corr(Gamma_RGDP(15:end-1),grRGDP_FRA(15:end-1))

yyaxis left
plot(year_WDI_(2:end-1),Gamma_RGDP(2:end-1),'r-');
ylabel('GDP Granular Residual')
yyaxis right
plot(year_WDI_(2:end-1),grRGDP_FRA(2:end-1),'b-');
ylabel('GDP Growth')
xlabel('')
saveas(gcf,'Figures/rho3eta1.000001lambda1.000001_psi3_sigma1.5/GDPGran_GDPGrowth_line.pdf')


yyaxis left
plot(year_WDI_(2:55),Gamma_TFP(12:end),'r-');
ylabel('TFP Granular Residual')
yyaxis right
plot(year_WDI_(2:55),grRGDP_FRA(2:55),'b-');
ylabel('GDP Growth')
xlabel('')
saveas(gcf,'Figures/rho3eta1.000001lambda1.000001_psi3_sigma1.5/TFPGran_GDPGrowth_line.pdf')


% Compute granular residual for France based on results of the Econometrica
% paper

load('MAT/External/granular_residual');

b = regstats(Gamma_TFP(42:58),agr_idio,'linear','yhat');
labels = cellstr(num2str(year_ficus));
figure(1)
plot(agr_idio,Gamma_TFP(42:58),'bo',agr_idio,b.yhat,'r-','LineWidth',2)
%text(agr_idio,Gamma_TFP(42:58),labels,'VerticalAlignment','bottom','HorizontalAlignment','right')
xlabel('French Granular Residual')
ylabel('TFP Granular Residual')
saveas(gcf,'Figures/rho3eta1.000001lambda1.000001_psi3_sigma1.5/TFPGran_vs_FrenchGran.pdf')
clear b
std(grRGDP_FRA(32:48))
std(agr_idio)
std(Gamma_TFP(42:58))
corr(Gamma_TFP(42:58),grRGDP_FRA(32:48))
corr(Gamma_TFP(42:58),agr_idio)

b = regstats(Gamma_RGDP(32:48),agr_idio,'linear','yhat');
figure(2)
plot(agr_idio,Gamma_RGDP(32:48),'bo',agr_idio,b.yhat,'r-','LineWidth',2)
%text(agr_idio,Gamma_RGDP(32:48),labels,'VerticalAlignment','bottom','HorizontalAlignment','right')
xlabel('French Granular Residual')
ylabel('GDP Granular Residual')
saveas(gcf,'Figures/rho3eta1.000001lambda1.000001_psi3_sigma1.5/GDPGran_vs_FrenchGran.pdf')
clear b
std(grRGDP_FRA(32:48))
std(agr_idio)
std(Gamma_RGDP(32:48))
corr(Gamma_RGDP(32:48),grRGDP_FRA(32:48))
corr(Gamma_RGDP(32:48),agr_idio)

%% Generate table of output using write table

% Data moments

sdRGDPgr_7514 = std(grRGDP_FRA(year_WDI_>=1975&year_WDI_<=2014))*100;
sdRGDPgr_9107 = std(grRGDP_FRA(year_WDI_>=1991&year_WDI_<=2007))*100;
stGRANres_9107 = sd_agidio*100;

% Foreign TFP moments

sddlYF_TFP_7514 = std(Overall_TFP(year_pwt_>=1975&year_pwt_<=2014))*100/2;
sddGF_TFP_7514 = std(Gamma_TFP(year_pwt_>=1975&year_pwt_<=2014))*100/2;
sdEF_TFP_7514 = std(Mean_TFP(year_pwt_>=1975&year_pwt_<=2014))*100/2;

sddlYF_TFP_9107 = std(Overall_TFP(year_pwt_>=1991&year_pwt_<=2007))*100/2;
sddGF_TFP_9107 = std(Gamma_TFP(year_pwt_>=1991&year_pwt_<=2007))*100/2;
sdEF_TFP_9107 = std(Mean_TFP(year_pwt_>=1991&year_pwt_<=2007))*100/2;

% Foreign GDP moments

sddlYF_RGDP_7514 = std(Overall_RGDP(year_WDI_>=1975&year_WDI_<=2014))*100;
sddGF_RGDP_7514 = std(Gamma_RGDP(year_WDI_>=1975&year_WDI_<=2014))*100;
sdEF_RGDP_7514 = std(Mean_RGDP(year_WDI_>=1975&year_WDI_<=2014))*100;

sddlYF_RGDP_9107 = std(Overall_RGDP(year_WDI_>=1991&year_WDI_<=2007))*100;
sddGF_RGDP_9107 = std(Gamma_RGDP(year_WDI_>=1991&year_WDI_<=2007))*100;
sdEF_RGDP_9107 = std(Mean_RGDP(year_WDI_>=1991&year_WDI_<=2007))*100;

% Collect terms for table

table_data = [sdRGDPgr_7514 NaN;sdRGDPgr_9107 stGRANres_9107];
table_tfp = [sddlYF_TFP_7514  sdEF_TFP_7514  sddGF_TFP_7514;sddlYF_TFP_9107  sdEF_TFP_9107 sddGF_TFP_9107];
table_rgdp = [sddlYF_RGDP_7514  sdEF_RGDP_7514 sddGF_RGDP_7514; sddlYF_RGDP_9107  sdEF_RGDP_9107 sddGF_RGDP_9107];

table_granres = table(table_data,table_tfp,table_rgdp);
table_granres.Properties.RowNames = {'1975--2014' '1991--2007'};

writetable(table_granres,'Tables/table4_granres.xlsx','WriteVariableNames',true,'WriteRowNames',true);
